/*
 * ULSF Sagittarius IDA Jones Panel
 * Copyright (C) 2002, Kevin Jones, all rights reserved.
 * Copyright (C) 2012, John Pritchard, all rights reserved.
 * 
 * This program is free software: you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see
 * <http://www.gnu.org/licenses/>.
 */
package sagittarius.panel;

import lxl.Map;

import java.io.IOException;
import java.io.FileWriter;
import java.io.PrintWriter;
import java.net.URL;
import java.net.MalformedURLException;

/**
 * Java port of Kevin's Online Panel Code.
 * @see http://aa.nps.edu/~jones/online_tools/panel2
 */
public class Panel {

    public final static float PI = (float)Math.PI;
    public final static float PI2INV = (0.5f / PI);

    /**
     * After the panel has been initialized, the "xb" and "yb" arrays
     * are coordinates between zero and one, from the trailing edge at
     * (1,0) to the leading edge at (0,0) in clockwise order --
     * subject to any non zero angle of attack input to the panel
     * constructor.
     * 
     * The "xm" and "ym" are midpoints for the panel between
     * coordinates n and n+1.
     * 
     * The leading edge begins with index "nle".
     * 
     * <h3>Dimensions</h3>
     * 
     * Note that the traditional cartesian (X,Y) notation employed
     * here maps to the traditional aircraft "X" and "Z" dimensions.
     */
    public final static class Surface {

        /**
         * Computational cache object serves as lookup key and value
         */
        public final static class Area
            extends java.lang.Object
            implements java.lang.Comparable<Area>
        {

            public final float scaleX, scaleY;

            public final float area;

            public final String string;

            public final int hashCode;


            public Area(float sx, float sy, float a){
                super();
                this.scaleX = sx;
                this.scaleY = sy;
                this.area = a;

                this.string = String.format("(%g %g) %g",sx,sy,a);
                /*
                 * See Float.hashCode
                 */
                {
                    final int hashX = Float.floatToIntBits(sx);
                    final int hashY = Float.floatToIntBits(sy);
                    this.hashCode = (hashX^hashY);
                }
            }

            public int compareTo(Area that){

                if (this.scaleX == that.scaleX){

                    if (this.scaleY == that.scaleY)

                        return 0;

                    else if (this.scaleY < that.scaleY)
                        return -1;
                    else
                        return 1;
                }
                else if (this.scaleX < that.scaleX)
                    return -1;
                else
                    return 1;
            }
            public int hashCode(){
                return this.hashCode;
            }
            public String toString(){
                return this.string;
            }
            public boolean equals(Object that){
                if (this == that)
                    return true;

                else if (that instanceof Area){
                    Area thatA = (Area)that;
                    return (this.scaleX == thatA.scaleX && this.scaleY == thatA.scaleY);
                }
                else
                    return false;
            }
        }

        /**
         * Computational cache object serves as lookup key and value
         */
        public final static class Centroid
            extends java.lang.Object
            implements java.lang.Comparable<Centroid>
        {

            public final float scaleX, scaleY;

            public final float x, y;

            public final String string;

            public final int hashCode;


            public Centroid(float sx, float sy, float cx, float cy){
                super();
                this.scaleX = sx;
                this.scaleY = sy;
                this.x = cx;
                this.y = cy;

                this.string = String.format("(%g %g) %g %g",sx,sy,cx,cy);
                /*
                 * See Float.hashCode
                 */
                {
                    final int hashX = Float.floatToIntBits(sx);
                    final int hashY = Float.floatToIntBits(sy);
                    this.hashCode = (hashX^hashY);
                }
            }

            public int compareTo(Centroid that){

                if (this.scaleX == that.scaleX){

                    if (this.scaleY == that.scaleY)

                        return 0;

                    else if (this.scaleY < that.scaleY)
                        return -1;
                    else
                        return 1;
                }
                else if (this.scaleX < that.scaleX)
                    return -1;
                else
                    return 1;
            }
            public int hashCode(){
                return this.hashCode;
            }
            public String toString(){
                return this.string;
            }
            public boolean equals(Object that){
                if (this == that)
                    return true;

                else if (that instanceof Centroid){
                    Centroid thatC = (Centroid)that;
                    return (this.scaleX == thatC.scaleX && this.scaleY == thatC.scaleY);
                }
                else
                    return false;
            }
        }

        public final int npoints, nle;
        public float xb[], yb[], xm[], ym[], dx[], dy[], cos_theta[], sin_theta[];
        public float cp[], cpmin, cpmax, surface;

        private final Map<Area,Area> area = new Map<Area,Area>();
        private final Map<Centroid,Centroid> centroid = new Map<Centroid,Centroid>();


        public Surface(Panel p){
            super();
            this.npoints = (p.inputs.npanels+1);
            this.nle = (p.inputs.npanels/2);
            this.xb = new float[p.npdim];
            this.yb = new float[p.npdim];
            this.xm = new float[p.npdim-1];
            this.ym = new float[p.npdim-1];
            this.dx = new float[p.npdim-1];
            this.dy = new float[p.npdim-1];
            this.cos_theta = new float[p.npdim];
            this.sin_theta = new float[p.npdim];
            this.cp = new float[p.npdim-1];
        }

        public void print(){
            PrintWriter out = null;
            try {
                out = new PrintWriter(new FileWriter("airfoil.fmt"));
                for (int idx = 0, len = this.npoints; idx < len; idx++){
                    out.printf("  %g  %g\n",this.xb[idx],this.yb[idx]);
                }
            }
            catch (IOException exc){
                exc.printStackTrace();
            }
            finally {
                if (null != out)
                    out.close();
            }
        }
        /**
         * @return Coordinate index of trailing edge
         */
        public int indexTE(){
            return 0;
        }
        /**
         * @return Coordinate index of leading edge
         */
        public int indexLE(){
            return this.nle;
        }
        public int coordCount(){
            return this.npoints;
        }
        public float coordX(int index, float scale){
            return (this.xb[index]*scale);
        }
        public float coordY(int index, float scale){
            return (this.yb[index]*scale);
        }
        /**
         * @return The negative area of the clockwise polygon.
         * @see java.lang.Math#abs(float)
         */
        public Area area(float scaleX, float scaleY){

            if (0.0f < scaleX && 0.0f < scaleY){

                Area area = this.area.get(new Area(scaleX,scaleY,0));

                if (null == area){
                    float A = 0.0f;

                    for (int cc = 0, zz = (this.npoints-2); cc < zz; cc++){
                        float x0 = this.xb[cc]*scaleX;
                        float x1 = this.xb[cc+1]*scaleX;

                        float y0 = this.yb[cc]*scaleY;
                        float y1 = this.yb[cc+1]*scaleY;

                        A += (x0*y1 - x1*y0);
                    }
                    A /= 2.0f;

                    area = new Area(scaleX,scaleY,A);

                    this.area.put(area,area);
                }
                return area;
            }
            else
                throw new IllegalArgumentException(String.format("(scaleX %g) (scaleY %g)",scaleX,scaleY));
        }
        /**
         * @return (X,Y)
         */
        public Centroid centroid(float scaleX, float scaleY){

            if (0.0f < scaleX && 0.0f < scaleY){

                Centroid centroid = this.centroid.get(new Centroid(scaleX,scaleY,0,0));

                if (null == centroid){

                    final float area6 = (6.0f * (this.area(scaleX,scaleY).area));

                    float cx = 0.0f, cy = 0.0f;

                    for (int cc = 0, zz = (this.npoints-2); cc < zz; cc++){
                        float x0 = this.xb[cc]*scaleX;
                        float x1 = this.xb[cc+1]*scaleX;

                        float y0 = this.yb[cc]*scaleY;
                        float y1 = this.yb[cc+1]*scaleY;

                        float f = ((x0*y1)-(x1*y0));

                        cx += ((x0+x1)*f);

                        cy += ((y0+y1)*f);
                    }

                    if (area6 == area6 && 0.0f != area6){

                        cx /= area6;

                        cy /= area6;
                    }
                    centroid = new Centroid(scaleX,scaleY,cx,cy);

                    this.centroid.put(centroid,centroid);
                }
                return centroid;
            }
            else
                throw new IllegalArgumentException(String.format("(scaleX %g) (scaleY %g)",scaleX,scaleY));
        }
    }
    /**
     * 
     */
    public final static class Field {

        public float x[][], y[][], pres[][], vel[][], u_vel[][], v_vel[][], xs[][], ys[][], ns[];

        public Field(Panel p){
            super();

            this.x = new float[p.npdim][p.kdim];
            this.y = new float[p.npdim][p.kdim];
            this.pres = new float[p.npdim][p.kdim];
            this.vel = new float[p.npdim][p.kdim];
            this.u_vel = new float[p.npdim][p.kdim];
            this.v_vel = new float[p.npdim][p.kdim];
            this.xs = new float[p.ldim][p.lsdim];
            this.ys = new float[p.ldim][p.lsdim];
            this.ns = new float[p.ldim];
        }
    }
    public final static class Singul {

        public float q[], gamma;

        public Singul(Panel p){
            super();

            this.q = new float[p.npdim];
        }
    }
    public final static class Coef {

        public float a[][];

        public Coef(Panel p){
            super();

            this.a = new float[p.npdim][p.npdim+1];
        }
    }
    public final static class Infcof {

        public float aan[][], bbn[][];

        public Infcof(Panel p){
            super();
            this.aan = new float[p.npdim][p.npdim];
            this.bbn = new float[p.npdim][p.npdim];
        }
    }
    public final static class Forces {

        public float cl, cd, cm;

        public Forces(Panel p){
            super();
        }

        public void print(){
            System.out.printf("%3.3f %3.3f %3.3f\n",this.cl,this.cd,this.cm);
        }
    }
    public final static class Consts {

        public final float alphar, cos_alpha, sin_alpha;
        public String title;
        public int icounter;
        public final float xmn, xmx, ymn, ymx;

        public Consts(Panel p){
            super();
            this.alphar = p.inputs.alphad * PI / 180.0f;
            this.cos_alpha = (float)Math.cos(this.alphar);
            this.sin_alpha = (float)Math.sin(this.alphar);
            this.xmn = p.inputs.x0 - p.inputs.width/2.0f;
            this.xmx = this.xmn + p.inputs.width;
            this.ymn = p.inputs.y0 - p.inputs.width/2.0f;
            this.ymx = this.ymn + p.inputs.width;
        }
    }
    public final static class Inputs {

        public final int naca;

        public final URL airfoil;

        public final int npanels;

        public final float alphad, xcp, ycp, width, x0, y0;

        public final boolean smoothing;

        public Inputs(int naca, boolean smoothing, int np, float ad, float xc, float yc, 
                      float w, float x0, float y0)
        {
            super();
            if (0 < naca && 0 < np && 0.0f < w){
                this.naca = naca;
                this.airfoil = null;
                this.smoothing= smoothing;
                this.npanels = ((np>>>1)<<1);//(ensure even)
                this.alphad = ad;
                this.xcp = xc;
                this.ycp = yc;
                this.width = w;
                this.x0 = x0;
                this.y0 = y0;
            }
            else
                throw new IllegalArgumentException();
        }
        public Inputs(URL airfoil, boolean smoothing, int np, float ad, float xc, float yc, 
                      float w, float x0, float y0)
        {
            super();
            if (null != airfoil && 0 < np && 0.0f < w){
                this.naca = 0;
                this.airfoil = airfoil;
                this.smoothing= smoothing;
                this.npanels = ((np>>>1)<<1);//(ensure even)
                this.alphad = ad;
                this.xcp = xc;
                this.ycp = yc;
                this.width = w;
                this.x0 = x0;
                this.y0 = y0;
            }
            else
                throw new IllegalArgumentException();
        }

        public boolean isNACA(){
            return (null == this.airfoil);
        }
        public boolean isDownload(){
            return (null != this.airfoil);
        }
    }


    public final static int npdim=1001, npdim1 = (npdim-1), npdim2 = (npdim1/2), kdim=31, ldim=41, lsdim=3000;

    public final Inputs inputs;

    public final Surface surface;

    public final Field field;

    public final Singul singul;

    public final Coef coef;

    public final Infcof infcof;

    public final Forces forces;

    public final Consts consts;

    public Airfoil airfoil;


    public Panel(int naca, boolean smoothing, int np, float ad, float xc, float yc, 
                 float w, float x0, float y0)
        throws IOException
    {
        this(new Inputs(naca, smoothing, Math.min(npdim2,np), ad, xc, yc, w, x0, y0));
    }
    public Panel(URL airfoil, boolean smoothing, int np, float ad, float xc, float yc, 
                 float w, float x0, float y0)
        throws IOException
    {
        this(new Inputs(airfoil, smoothing, Math.min(npdim2,np), ad, xc, yc, w, x0, y0));
    }
    public Panel(String airfoil, boolean smoothing, int np, float ad, float xc, float yc, 
                 float w, float x0, float y0)
        throws IOException
    {
        this(new Inputs(UIUC.For(airfoil), smoothing, Math.min(npdim2,np), ad, xc, yc, w, x0, y0));
    }
    public Panel(String airfoil)
        throws IOException
    {
        this(airfoil,false,npdim2,0,0,0,100,0,0);
    }
    /**
     * Load coordinates
     */
    public Panel(Inputs inputs)
        throws IOException
    {
        super();
        if (null != inputs){
            this.inputs = inputs;
            this.surface = new Surface(this);
            this.field = new Field(this);
            this.singul = new Singul(this);
            this.coef = new Coef(this);
            this.infcof = new Infcof(this);
            this.forces = new Forces(this);
            this.consts = new Consts(this);

            this.loadCoordinates();
        }
        else
            throw new IllegalArgumentException();
    }


    public boolean isNACA(){
        return this.inputs.isNACA();
    }
    public boolean isDownload(){
        return this.inputs.isDownload();
    }
    /**
     * Compute Coefficients and CP
     */
    public void init()
    { 
        this.loadCoefficients();
        this.gaussianElimination();
        this.loadCp();
    }
    private void loadCoordinates()
        throws IOException
    {
        if (this.inputs.isNACA())
            this.airfoil = new NACA(this);
        else
            this.airfoil = new UIUC(this);

        this.surface.surface = 0.0f;
        float[] xm = this.surface.xm;
        float[] ym = this.surface.ym;
        float[] xb = this.surface.xb;
        float[] yb = this.surface.yb;
        float[] dx = this.surface.dx;
        float[] dy = this.surface.dy;
        float[] sin_theta = this.surface.sin_theta;
        float[] cos_theta = this.surface.cos_theta;
        float ds;

        for (int n = 0, z = this.inputs.npanels; n < z; n++){
            /*
             * Panel mid points
             */
            xm[n] = (xb[n+1] + xb[n])/2.0f;
            ym[n] = (yb[n+1] + yb[n])/2.0f;
            /*
             * Arc length
             */
            dx[n] = xb[n+1] - xb[n];
            dy[n] = yb[n+1] - yb[n];

            ds = (float)Math.sqrt( Math.pow(dx[n],2) + Math.pow(dy[n],2) );

            this.surface.surface += ds;
            /*
             * Slope
             */
            sin_theta[n] = dy[n] / ds;
            cos_theta[n] = dx[n] / ds;
        }
    }
    private void loadCoefficients(){
        final int z = this.surface.npoints, y = (z-1);

        float[][] a = this.coef.a;
        for (int n = 0; n < z; n++){
            a[y][n] = 0.0f;
        }

        float flog, ftan, dxj, dxjp, dyj, dyjp, ctimtj, stimtj;

        float[][] aan = this.infcof.aan, bbn = this.infcof.bbn;

        float[] xm = this.surface.xm;
        float[] ym = this.surface.ym;
        float[] xb = this.surface.xb;
        float[] yb = this.surface.yb;
        float[] cos_theta = this.surface.cos_theta;
        float[] sin_theta = this.surface.sin_theta;
        float cos_alpha = this.consts.cos_alpha;
        float sin_alpha = this.consts.sin_alpha;

        for (int i = 0, j; i < z; i++){
            a[i][y] = 0.0f;

            for (j = 0; j < z; j++){

                if (i == j){
                    flog = 0.0f;
                    ftan = 0.0f;
                }
                else {
                    dxj  = xm[i] - xb[j];
                    dxjp = xm[i] - xb[j+1];
                    dyj  = ym[i] - yb[j];
                    dyjp = ym[i] - yb[j+1];
                    flog = (float)Math.log( ( dxjp*dxjp + dyjp*dyjp )
                                     / ( dxj*dxj  + dyj*dyj  )) / 2;
                    ftan = (float)Math.atan2( dyjp * dxj - dxjp * dyj,
                                       dxjp * dxj + dyjp * dyj );
                }

                ctimtj = (cos_theta[i] * cos_theta[j]
                          + sin_theta[i] * sin_theta[j]);
                stimtj = (sin_theta[i] * cos_theta[j]
                          - cos_theta[i] * sin_theta[j]);
                
                a[i][j] = PI2INV * ( ftan * ctimtj + flog * stimtj );
                bbn[i][j] = PI2INV * ( flog * ctimtj - ftan * stimtj );

                a[i][y] = a[i][y] + bbn[i][j];

                if (( 1 == i) || ( y == i)){
                    a[y][j] = a[y][j] - bbn[i][j];
                    a[y][y] = a[y][y] + a[i][j];
                }

                aan[i][j] = a[i][j];
            }

            a[i][z] = (sin_theta[i] * cos_alpha
                       - cos_theta[i] * sin_alpha);
        }

        a[y][z] = (-( cos_theta[0] + cos_theta[y])
                   * cos_alpha
                   -( sin_theta[0] + sin_theta[y])
                   * sin_alpha);
    }
    private void gaussianElimination(){
        final int z = this.surface.npoints, y = (z-1);

        float[][] a = this.coef.a;
        float r;

        for (int h = 0, i = 1, j, k; i < z; h++, i++){

            for (j = i; j < z; j++){

                r = a[j][h] / a[h][h];

                for (k = i; k <= z; k++){

                    a[j][k] = a[j][k] - r * a[h][k];
                }
            }
        }

        a[y][z] = a[y][z] / a[y][y];

        for (int l = 1, i, j; l < y; l++){

            i = y + 1 - l;

            for ( j = i+1; j < y; j++){

                a[i][z] = a[i][z] - a[i][j] * a[j][z];
            }

            a[i][z] = a[i][z] / a[i][i];
        }
    }
    private void loadCp(){
        final int x = this.inputs.npanels, w = (x-1), z = this.surface.npoints, y = (z-1);

        float[] q = this.singul.q;
        final float[][] a = this.coef.a;

        for (int i = 0; i < x; i++){

            q[i] = a[i][z];
        }

        final float gamma = this.singul.gamma = a[y][z];
        final float[] sin_theta = this.surface.sin_theta;
        final float[] cos_theta = this.surface.cos_theta;
        final float[] xm = this.surface.xm;
        final float[] ym = this.surface.ym;
        final float[] dx = this.surface.dx;
        final float[] dy = this.surface.dy;
        final float cos_alpha = this.consts.cos_alpha;
        final float sin_alpha = this.consts.sin_alpha;
        final float[][] aan = this.infcof.aan, bbn = this.infcof.bbn;
        final float xcp = this.inputs.xcp;
        final float ycp = this.inputs.ycp;
        float[] cp = this.surface.cp;

        this.surface.cpmin = 1.0f;
        this.surface.cpmax = -1.0f;
        this.forces.cm = 0.0f;

        float cfx = 0.0f, cfy = 0.0f;
        float vtfree, vtang, xmid, ymid;

        for (int i = 0, j; i < x; i++){

            vtfree = cos_alpha * cos_theta[i] + sin_alpha * sin_theta[i];
            vtang = vtfree;

            for (j = 0; j < x; j++){

                vtang = vtang - bbn[i][j] * q[j] + gamma * aan[i][j];
            }
            cp[i] = 1.0f - vtang*vtang;

            if ( cp[i] < this.surface.cpmin )
                this.surface.cpmin = cp[i];

            if ( cp[i] > this.surface.cpmax )
                this.surface.cpmax = cp[i];

            xmid = xm[i] - xcp;
            ymid = ym[i] - ycp;
            cfx = cfx + cp[i] * dy[i];
            cfy = cfy - cp[i] * dx[i];

            this.forces.cm = this.forces.cm + cp[i] * ( dx[i] * xmid + dy[i] * ymid );
       }

        this.forces.cl = cfy * cos_alpha - cfx * sin_alpha;
        this.forces.cd = cfx * cos_alpha + cfy * sin_alpha;
    }
    public void loadGrid(){
        final int kmax = this.kdim;

        final int x = this.inputs.npanels, w = (x-1), z = this.surface.npoints, y = (z-1);
        final int nle = this.surface.nle;

        final float[] xb = this.surface.xb;
        final float[] yb = this.surface.yb;

        float[][] gridX = this.field.x;
        float[][] gridY = this.field.y;

        final float xmn = this.consts.xmn, xmx = this.consts.xmx, ymn = this.consts.ymn, ymx = this.consts.ymx;
        final float rad1 = (float)Math.sqrt(Math.pow((xmn-0.5f),2)+(ymn*ymn));
        final float rad2 = (float)Math.sqrt(Math.pow((xmx+0.5f),2)+(ymn*ymn));
        final float rad3 = (float)Math.sqrt(Math.pow((xmx+0.5f),2)+(ymx*ymx));
        final float rad4 = (float)Math.sqrt(Math.pow((xmn-0.5f),2)+(ymx*ymx));
        final float radius = Math.max(0.75f,Math.max(Math.max(rad1,rad2),Math.max(rad3,rad4)));
        final float ds_wall = 1.0e-4f;

        final float xte = ( xb[0] + xb[y]) / 2.0f;
        final float yte = ( yb[0] + yb[y]) / 2.0f;

        float dxb = xte - ( xb[1] + xb[y-1]) / 2.0f;
        float dyb = yte - ( yb[1] + yb[y-1]) / 2.0f;

        float slope = (dyb / dxb), theta, theta_te, theta_le, b_int, aq, bq, cq, radical, xo, yo, ds;

        if ( slope < 1.0e-6f )
            theta_te = 0.0f;
        else {
            b_int = yte - slope * xte;

            aq = 1.0f + slope*slope;
            bq = 2.0f * slope * b_int - 1.0f;
            cq = b_int*b_int - radius*radius + 0.25f;

            radical = bq*bq - 4.0f * aq * cq;

            if ( radical < 0.0f )
                throw new IllegalStateException("Negative radical");
            else {
                xo = ( -bq + (float)Math.sqrt( radical )) / ( 2.0f * aq );
                yo = b_int + slope * xo;

                theta_te = (float)Math.atan2(yo,xo-0.5f);
            }
        }

        dxb = xb[nle+1] - xb[nle-1];
        dyb = yb[nle+1] - yb[nle-1];

        slope = -dxb / dyb;

        if (( Math.abs(dxb) < 1.0e-6f ) ||
            ( Math.abs(slope) < 1.0e-6f ))
        {
            theta_le = -PI;
        }
        else {
            b_int = yb[nle] - slope * xb[nle];

            aq = 1.0f + slope*slope;
            bq = 2.0f * slope * b_int - 1.0f;
            cq = b_int*b_int - radius*radius + 0.25f;

            radical = bq*bq - 4.0f * aq * cq;
            if ( radical < 0.0f )
                throw new IllegalStateException("Negative radical");
            else {

                xo = ( -bq - (float)Math.sqrt( radical )) / ( 2.0f * aq );
                yo = b_int + slope * xo;

                theta_le = (float)Math.atan2(yo,xo-0.5f);
                if ( theta_le > 0.0 ) 
                    theta_le = theta_le - 2.0f*PI;
            }
        }

        final float dtheta_lower = ( theta_le - theta_te ) / (nle-1);
        final float dtheta_upper = ( theta_te - theta_le-2.0f*PI)/(this.inputs.npanels-nle+1);

        float[] ct = new float[10], c = new float[10];
        int[] ipt = new int[10], ict = new int[10];
        float dgx, dgy, s;

        for (int j = 0, nle1 = (nle-1); j < nle1; j++){

            gridX[j][0] = xb[j];
            gridY[j][0] = yb[j];

            theta = theta_te + dtheta_lower * ( j-1 );
            gridX[j][kmax-1] = 0.5f + radius * (float)Math.cos(theta);
            gridY[j][kmax-1] = radius * (float)Math.sin(theta);

            dgx = gridX[j][kmax-1] - gridX[j][0];
            dgy = gridY[j][kmax-1] - gridY[j][0];
            ds = (float)Math.sqrt( dgx*dgx + dgy*dgy );

            ipt[0] = 1;
            ict[0] = 0;
            ct[0] = 0.0f;
            ipt[1] = 2;
            ict[1] = 0;
            ct[1] = ds_wall;
            ipt[2] = kmax;
            ict[2] = 0;
            ct[2] = ds;

            if (Airfoil.Polynomial(ct,ict,ipt,3,c))
                throw new IllegalStateException("Singular s-redistribution");
            else {
                for (int k = 1, kmax1 = (kmax-1); k < kmax1; k++){

                    s = c[0] + c[1]*k + (float)Math.pow(c[2]*k,2);
                    gridX[j][k] = s/ds * dgx + gridX[j][0];
                    gridY[j][k] = s/ds * dgy + gridY[j][0];
                }
            }
        }

        for (int j = nle; j < x; j++){

            gridX[j][0] = xb[j];
            gridY[j][0] = yb[j];

            theta = theta_le + dtheta_upper * ( j-nle );
            gridX[j][kmax-1] = 0.5f + radius * (float)Math.cos(theta);
            gridY[j][kmax-1] = radius * (float)Math.sin(theta);
            
            dgx = gridX[j][kmax-1] - gridX[j][0];
            dgy = gridY[j][kmax-1] - gridY[j][0];
            ds = (float)Math.sqrt( dgx*dgx + dgy*dgy );

            ipt[0] = 1;
            ict[0] = 0;
            ct[0] = 0.0f;
            ipt[1] = 2;
            ict[1] = 0;
            ct[1] = ds_wall;
            ipt[2] = kmax;
            ict[2] = 0;
            ct[2] = ds;

            if (Airfoil.Polynomial(ct,ict,ipt,3,c))
                throw new IllegalStateException("Singular s-redistribution");
            else {

                for (int k = 1, kmax1 = (kmax-1); k < kmax1; k++){

                    s = c[0] + c[1]*k + c[2]*k*k;
                    gridX[j][k] = s/ds * dgx + gridX[j][0];
                    gridY[j][k] = s/ds * dgy + gridY[j][0];
                }
            }
        }

        for (int k = 1; k < kmax; k++){

            gridX[y][k] = gridX[0][k];
            gridY[y][k] = gridY[0][k];
        }
    }
    public void loadFieldVars(){
        final int jmax = this.surface.npoints, jmax1 = (jmax-1);
        final int kmax = this.kdim;
        final float x[][] = this.field.x, y[][] = this.field.y;
        float pres[][] = this.field.pres;
        float vel[][] = this.field.vel, u_vel[][] = this.field.u_vel, v_vel[][] = this.field.v_vel;

        float vsq;
        Airfoil.Velocity uv = new Airfoil.Velocity();

        for (int j = 0, k; j < jmax1; j++){

            for (k = 1; k < kmax; k++){

                Airfoil.Velocity(this,x[j][k],y[j][k],uv);
                u_vel[j][k] = uv.u;
                v_vel[j][k] = uv.v;
                vsq = uv.u*uv.u + uv.v*uv.v;
                vel[j][k] = (float)Math.sqrt(vsq);
                pres[j][k] = 1.0f - vsq;
            }
        }

        for (int k = 1; k < kmax; k++){

            u_vel[jmax1][k] = u_vel[0][k];
            v_vel[jmax1][k] = v_vel[0][k];
            vel[jmax1][k]   = vel[0][k];
            pres[jmax1][k]  = pres[0][k];
        }
    }
    public void loadStreamLines(){
        final float xmx = this.consts.xmx;
        final float xmn = this.consts.xmn;
        final float ymx = this.consts.ymx;
        final float ymn = this.consts.ymn;

        float xs[][] = this.field.xs, ys[][] = this.field.ys, ns[] = this.field.ns;

        /*
         * Integration time-step.
         */
        final float dt = (xmx-xmn)/200;

        final float xstart = Math.min(-0.01f,xmn);

        final float alphar2 = (this.forces.cl/PI+this.consts.alphar)/2.0f;

        final float slope = (float)Math.tan(alphar2);
        final float b_int = this.inputs.y0 - slope * (this.inputs.x0+0.25f)/2;

        final float y_cent = xstart * slope + b_int, y_inc;

        if ( y_cent < this.inputs.y0 )
            y_inc = 2.0f * ( ymx - y_cent ) / this.ldim;
        else
            y_inc = 2.0f * ( y_cent - ymn ) / this.ldim;

        float dx1, dy1, dx2, dy2, dx3, dy3, dx4, dy4;

        Airfoil.Velocity uv = new Airfoil.Velocity();

        for (int l = 0, ldim = this.ldim, ldim2 = (ldim/2), m; l < ldim; l++){

            m = 0;
            xs[l][m] = xstart;
            ys[l][m] = y_cent + y_inc * (l-ldim2-1);
            /*
             * Use RK-4th-Order until we leave the domain or we hit mdim points.
             */
            do {
                Airfoil.Velocity(this,xs[l][m],ys[l][m],uv);
                dx1 = dt * uv.u;
                dy1 = dt * uv.v;

                Airfoil.Velocity(this,xs[l][m]+dx1/2,ys[l][m]+dy1/2,uv);
                dx2 = dt * uv.u;
                dy2 = dt * uv.v;

                Airfoil.Velocity(this,xs[l][m]+dx2/2,ys[l][m]+dy2/2,uv);
                dx3 = dt * uv.u;
                dy3 = dt * uv.v;

                Airfoil.Velocity(this,xs[l][m]+dx3,ys[l][m]+dy3,uv);
                dx4 = dt * uv.u;
                dy4 = dt * uv.v;

                xs[l][m+1] = xs[l][m] + ( dx1 + dx2*2 + dx3*2 + dx4) / 6;
                ys[l][m+1] = ys[l][m] + ( dy1 + dy2*2 + dy3*2 + dy4) / 6;

                m += 1;
            }
            while (( xs[l][m] < xmx ) && ( m < this.lsdim ));
            /*
             * Set the streamline length for drawing.
             */
            ns[l] = m;
        }
    }
    public void printForces(){
        this.forces.print();
    }
    public void printFoil(){
        this.surface.print();
    }
    public void printFieldPres(){
    }
    public void printGrid(){
    }
    public void printVectors(){
    }
    public void printStreamLines(){
    }
    public void printCp(){
    }

    public final static float Z(float v){
        if (1.0e-3f > Math.abs(v))
            return 0.0f;
        else
            return v;
    }
}
